home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
Software Vault: The Gold Collection
/
Software Vault - The Gold Collection (American Databankers) (1993).ISO
/
cdr53
/
acmalg01.zip
/
ACM615.FOR
< prev
next >
Wrap
Text File
|
1993-01-01
|
55KB
|
1,819 lines
C ALGORITHM 615, COLLECTED ALGORITHMS FROM ACM.
C ALGORITHM APPEARED IN ACM-TRANS. MATH. SOFTWARE, VOL.10, NO. 2,
C JUN., 1984, P. 202-206.
C PROGRAM SUBL1 (INPUT,OUTPUT,TAPE5=INPUT,TAPE6=OUTPUT)
C****************************************************************
C *THIS CODE SOLVES THE L1 BEST SUBSET PROBLEM*
C *SUBMITTED TO ACM TRANS.SPRING 1982*
C *LAST REVISED MAR. 12, 1984*
C****************************************************************
DIMENSION X(300,20), Y(300)
DIMENSION ZL(20), BVAL(210), IDEX(210), ISTAT(20)
C
COMMON /BLOCK/
- ACU,BIG,BSAVE(210),LU(20,20),PIE(210),SAD,SADT,SAVTOT(210),
-SIG(6000),SIGMA(300),ZLOW,ZSAVE(20),IBASE(20),IK,IK1,INCOL(20),
-INDEX(20),INPROB(210),IPAR(20),LABEL(210),LEVEL,DIRECT,INTL
C
INTEGER I,IBASE,IDEX,IFAULT,IK,IK1,IMT,INCOL,INDEX,INPROB,IPAR
INTEGER ISTAT,ITER,J,JMIN,K,L,LABEL,LEVEL,M,MININ,MMAX
INTEGER N,NAME,NMAX,NPROB
C
DOUBLE PRECISION ACU,BIG,BSAVE,BVAL,LU,PIE,POPT,SAD,SADT,SAVTOT
C++ REAL ACU,BIG,BSAVE,BVAL,LU,PIE,POPT,SAD,SADT,SAVTOT
DOUBLE PRECISION SIG,SIGMA,X,Y,ZL,ZLOW,ZSAVE
C++ REAL SIG,SIGMA,X,Y,ZL,ZLOW,ZSAVE
C
C
LOGICAL DIRECT,INTL
INTEGER UNITGO,UNITIN
C
C-----UNITIN IS THE INPUT UNIT NUMBER
C-----UNITGO IS THE OUTPUT UNIT NUMBER
C
DIMENSION IMT(72), NAME(40)
C
C-----IMT IS USED TO STORE THE INPUT FORMAT
C-----NAME SAVES THE PROBLEM NAME
C
UNITIN=5
UNITGO=6
NMAX=300
MMAX=30
NPROB=0
C
C-----K STORES THE NUMBER OF PARAMETERS IN THE FULL MODEL
C
10 READ (UNITIN,50,END=45) K,NAME
WRITE(UNITGO,50)K,NAME
IF (K.EQ.0) STOP
WRITE (UNITGO,60) NAME
C
C-----N STORES THE NUMBER OF OBSERVATIONS
C-----MININ STORES THE MINIMUM NUMBER OF PARAMETERS CONSIDERED
C-----POPT DEVIATION FROM OPTIMALITY ALLOWED
C
READ (UNITIN,130) N,MININ,POPT
51 FORMAT (4X,22HNUMBER OF OBSERVATION=,I5/4X,21HNUMBER OF PARAMETERS
1=,I5/4X,40HMINIMUM NUMBER OF PARAMETERS CONSIDERED=,I5/4X,45HPERCE
2NTAGE DEVIATION FROM OPTIMALITY ALLOWED=,F6.2/8X,2H**,
3 23HBEST SUBSET LAV PROGRAM)
WRITE (UNITGO,51) N,K,MININ,POPT
READ (UNITIN,150) IMT
DO 20 I=1,N
READ (UNITIN,IMT) Y(I), (X(I,J),J=1,K)
C
C WRITE (UNITGO,IMT) Y(I), (X(I,J),J=1,K)
C
20 CONTINUE
C
READ (UNITIN,70) (ISTAT(I),I=1,K)
C
ITER=0
C
C******CALL THE ROUTINE TO FIND THE BEST SUBSETS
C
C
CALL KBEST (X,Y,K,N,ITER,IFAULT,POPT,MININ,NMAX,MMAX,BVAL,IDEX,
1ISTAT,ZL)
C
WRITE (UNITGO,80) IFAULT
C
C WRITE THE FINAL BEST SUBSET SOLUTION
C
J=K*(K+1)/2
JMIN=(MININ-1)*(MININ)/2
J=J-JMIN
DO 40 I=MININ,K
WRITE (UNITGO,110) I
WRITE (UNITGO,100) ZL(I)
DO 30 L=1,I
WRITE (UNITGO,90) IDEX(J),BVAL(J)
J=J-1
30 CONTINUE
40 CONTINUE
C
NPROB=NPROB+1
WRITE (UNITGO,160)ITER
GO TO 10
C
45 STOP
C
50 FORMAT (I5,3X,40A1)
60 FORMAT (3X,5H ,14HPROBLEM TITLE ,2X,40A1)
70 FORMAT (10I2)
80 FORMAT (8H IFAULT=,I3)
90 FORMAT (3X/6H BETA(,I3,1H),F15.3)
100 FORMAT (4X,14(1H*),23HSUM OF ABSOLUTE VALUES=,F15.3)
110 FORMAT (//3X,36HBEST RESULTS FOR LAV SUBSET OF SIZE=,I3)
130 FORMAT (I5,I2,F6.2)
150 FORMAT (72A1)
160 FORMAT (//5X,16HITERATION COUNT=,I7)
C
END
SUBROUTINE KBEST (X,Y,M,N,ITER,IFAULT,POPT,MININ,NMAX,MMAX,BVAL
1,IDEX,ISTAT,ZL)
C
C***********************************************************************
C
C
C THE PURPOSE OF THIS PROGRAM IS TO DETERMINE THE BEST SUBSET OF
C PARAMETERS TO FIT A LINEAR REGRESSION UNDER AN LEAST ABSOLUTE
C VALUE CRITERION. THIS PROGRAM UTILIZES THE SIMPLEX METHOD OF
C LINEAR PROGRAMMING WITHIN A BRANCH-AND-BOUND ALGORITHM TO
C SOLVE THE BEST SUBSET PROBLEM.
C
C
C THE ALGORITHM IS BASED ON THE PUBLICATION:
C ARMSTRONG,R.D. AND M.T. KUNG "AN ALGORITHM TO SELECT THE BEST
C SUBSET FOR A LEAST ABSOLUTE VALUE REGRESSION PROBLEM,"
C OPTIMIZATION IN STATISTICS, TIMS STUDIES OF THE MANAGEMENT
C SCIENCES.
C
C FORMAL PARAMETERS
C
C X REAL ARRAY INPUT: VALUES OF INDEPENDENT VARIABLES
C (NMAX,MMAX) SUCH THAT EACH ROW CORREPSONDS TO
C AN OBSERVATION
C
C Y REAL ARRAY INPUT: VALUES OF THE DEPENDENT VARIABLES
C (NMAX)
C
C M INTEGER INPUT: NUMBER OF DEPENDENT VARIABLES
C
C N INTEGER INPUT: NUMBER OF OBSERVATIONS
C
C ITER INTEGER OUTPUT: NUMBER OF ITERATIONS
C
C IFAULT INTEGER OUTPUT: FAILURE INDICATOR
C =0 NORMAL TERMINATION
C =1 OBSERVATION MATRIX DOES NOT HAVE
C FULL ROW RANK (RANK M)
C =2 PROBLEM SIZE OUT OF RANGE
C =3 NO PIVOT ELEMENT FOUND IMPLYING
C NEAR SINGULAR BASIS
C
C POPT REAL INPUT: PERCENTAGE DEVIATION FROM
C OPTIMALITY ALLOWED
C
C MININ INTEGER INPUT: MINIMUM NUMBER OF PARAMETERS IN THE
C MODEL. BEST SUBSET OF SIZE MININ T
C M IS OBTAINED.
C
C NMAX INTEGER INPUT: DIMENSION OF ROWS IN X (ALSO Y)
C
C MMAX INTEGER INPUT: DIMENSION OF COLUMNS IN X
C
C BVAL REAL ARRAY OUTPUT: ARRAY OF OPTIMAL BETA VALUES FOR
C EACH SUBSET. THE BETA VALUES FOR
C THE SUBSET OF SIZE M ARE STORED
C IN POSITIONS BVAL(1),BVAL(2),...,
C BVAL(M), FOR THE SUBSET OF SIZE
C M-1 THE VALUES ARE STORED IN
C POSITIONS BVAL(M+1),BVAL(M+2),...,
C BVAL(2M-1). IN GENERAL, THE BETA
C VALUES FOR THE OPTIMAL SUBSET OF
C SIZE K ARE STORED IN POSITIONS
C BVAL(L),...,BVAL(L-K+1) WHERE
C L=(M*(M+1)-K*(K+1))/2 + 1
C
C IDEX INTEGER ARRAY OUTPUT: BETA INDEX SET FOR THE OPTIMAL
C (((MMAX+1)*MMAX)/2) SUBSET. THIS ARRAY IS A PARALLEL
C ARRAY FOR BVAL; I.E., IF BVAL(J)=2.
C AND IDEX(J)=7 THEN BETA(7)=2.7 IN
C THE ASSOCIATED OPTIMAL SUBSET.
C
C ISTAT INTEGER ARRAY INPUT: PARAMETER STATUS ARRAY.
C (MMAX)
C
C 1 IF BETA(J) IS REQUIRED
C IN EVERY MODEL
C ISTAT(J)=
C
C 0 OTHERWISE
C
C ZL REAL ARRAY OUTPUT: BEST OBJECTIVE VALUE FOR EACH SUBSE
C (MMAX)
C ZL(J) GIVES THE BEST OBJECTIVE VALU
C FOR THE SUBSET WITH M-J+1 PARAMETER
C***********************************************************************
C
C IMPLEMENTATION NOTES:
C
C 1. THE ROUTINE USES TWO MACHINE DEPENDENT VALUES ACU AND BIG.
C (THESE ARE SET AT THE BEGINNING OF THIS SUBROUTINE)
C A) ACU IS USED TO TEST FOR "ZERO". ACU SHOULD BE SET TO
C APPROXIMATELY 100 * THE RELATIVE MACHINE ACCURACY OF THE
C SYSTEM IN USE.
C B) BIG IS USED TO INITIALIZE THE ARRAY ZL (DESCRIBED ABOVE).
C BIG SHOULD BE SET TO THE LARGEST FLOATING POINT VALUE
C ASSIGNABLE.
C
C 2. BOTH SINGLE AND DOUBLE PRECISION VERSIONS ARE SUPPLIED. THIS
C VERSION IS IN DOUBLE PRECISION. TO GET A SINGLE PRECISION CO
C ALL STATEMENTS HAVING A "C++" IN COLUMNS 1-3 SHOULD BE MODIF
C BY DELETING THE "C++". THE CORRESPONDING DOUBLE PRECISION
C STATEMENTS SHOULD BE EITHER DELETED OR COMMENTED OUT.
C
C 3. ARRAY DIMENSIONS
C THE CODE IS CURRENTLY DIMENSIONED TO SOLVE PROBLEMS WITH UP
C 20 PARAMETERS AND 300 OBSERVATIONS. THE DIMENSION SIZES ARE
C DETERMINED AS FOLLOWS
C 20 MMAX
C 300 NMAX
C 210 (NMAX+1)*MMAX/2
C 6000 NMAX*MMAX
C 10 MAXIMUM VALUE OF IK .
C
C IF THE DIMENSION SIZES ARE CHANGED THE COMMON AND DIMENSION
C STATEMENTS WILL NEED TO BE MODIFIED IN EACH SUBROUTINE.
C
C
C***********************************************************************
C
C SUBROUTINES:
C
C CALBET - BACK-SOLVES A SYSYTEM OF EQUATIONS
C
C CALCPI - FORWARD-SOLVES A SYSYTEM OF EQUATIONS
C
C KBEST - THE DRIVER
C
C L1NORM - SOLVES THE INITIAL REGRESSION PROBLEM WITH ALL
C PARAMETERS INCLUDED IN THE MODEL
C
C PHASE2 - SOLVES THE CURRENT REGRESSION PROBLEM WITH A PRIMAL
C ALGORITHM
C
C SETUP - DETERMINES THE FORM OF THE CURRENT PROBLEM BY
C CHOOSING THE PARAMETER TO LEAVE BASED ON A PENALTY
C
C UPDATE - UPDATES LU DECOMPOSITION MATRIX
C
C
C***********************************************************************
C
C
C DESCRIPTION OF VARIABLES:
C
C IK:LENGTH OF CANDIDATE LIST
C BETA:ESTIMATES TO THE CURRENT SUBPROBLEM
C SAD:THE MINIMUM TOTAL ABSOLUTE DEVIATION
C IBASE:THE INDEX ARRAY OF COLUMNS OF THE BASIS
C
C N,M,NMAX,MMAX,X,Y, ARE NOT CHANGED IN THE SUBROUTINE;
C SAD IS UPDATED AT EACH ITERATION
C
C LU:THE LU DECOMPOSITION OF THE CURRENT BASIS
C INDEX:THE INDEX ARRAY OF ROWS OF LU
C TOT:THE CURRENT RHS OF THE DUAL PROBLEM
C SIGMA:INDICATOR ARRAY TO SPECIFY WHETHER A NONBASIC DUAL VARIABLE
C IS AT UPPER OR LOWER BOUND(+1 IMPLIES UPPER;-1 IMPLIES
C LOWER)
C INEXT:A LOCAL ARRAY FOR THE SORT ROUTINE
C RHS: A LOCAL ARRAY FOR THE CALBET ROUTINE
C
C IPAR(I)=-K IF THE K-TH PARAMETER (BETA) IS OUT OF MODEL AT LEVEL I
C = K IF THE K-TH PARAMETER IS IN
C
C LABEL(I) SAVES THE INDICES OF THE BASIC OBSERVATIONS IN THE
C PREDECESSOR PATH
C ZSAVE(I) SAVES THE OBJECTIVE VALUES ON THE PREDECESSOR PATH
C
C**********************************************************************
C
DIMENSION X(300,20), Y(300), ZL(20), BVAL(210), IDEX(210), ISTAT(2
10)
C
DIMENSION BETA(20), TOT(20), PI(20), REPP(20)
C
COMMON /BLOCK/
- ACU,BIG,BSAVE(210),LU(20,20),PIE(210),SAD,SADT,SAVTOT(210),
-SIG(6000),SIGMA(300),ZLOW,ZSAVE(20),IBASE(20),IK,IK1,INCOL(20),
-INDEX(20),INPROB(210),IPAR(20),LABEL(210),LEVEL,DIRECT,INTL
C
INTEGER I,IBASE,IDEX,IFAULT,IK,IK1,INCOL,INDEX,INPROB,IPAR
INTEGER ISTAT,ITER,J,JINM,K2,K3,KKK1,L,LABEL,LEVEL,M,MININ
INTEGER MMAX,N,NFREE,NMAX,NUMIN
C
DOUBLE PRECISION ACU,BETA,BIG,BSAVE,BVAL,LU,PI,PIE,POPT,POPT1
C++ REAL ACU,BETA,BIG,BSAVE,BVAL,LU,PI,PIE,POPT,POPT1
DOUBLE PRECISION REPP,SAD,SADT,SAVTOT,SIG,SIGMA,TOT,X,Y
C++ REAL REPP,SAD,SADT,SAVTOT,SIG,SIGMA,TOT,X,Y
DOUBLE PRECISION ZERO,ZL,ZLOW,ZSAVE
C++ REAL ZERO,ZL,ZLOW,ZSAVE
C
LOGICAL DIRECT,INTL
C
DATA ZERO/0.0D0/
C++ DATA ZERO/0.0E0/
C
C ASSIGN VALUES TO MACHINE DEPENDENT CONSTANTS
C SEE - FOX, P.A., A.D. HALL AND N.L. SCHRYER, "FRAMEWORK FOR
C A PORTABLE LIBRARY: ALGORITHM 528," ACM TRANSACTIONS
C ON MATHEMATICAL SOFTWARE, VOL 4, NO 2, JUNE 1978.
C
ACU = D1MACH(1)
ACU = 1000.*ACU
BIG = D1MACH(2)
C++ ACU = R1MACH(1)
C++ ACU = 100.*ACU
C++ BIG = R1MACH(2)
C
C
IFAULT=0
IK=5
IK1=IK-1
POPT1=(100.-POPT)/100.
C
C TEST PROBLEM SIZE
C
IF (N.LE.NMAX.AND.M.LE.MMAX) GO TO 10
IFAULT=2
RETURN
10 CONTINUE
DO 20 I=1,M
INDEX(I)=I
TOT(I)=ZERO
20 INCOL(I)=I
CALL L1NORM (X,Y,N,M,IFAULT,ITER,BETA,PI,TOT)
IF (IFAULT.GT.0) RETURN
C
C INITIALIZATION
C
C SAVE THE INITIAL SOLUTION
C
JINM=0
DO 30 I=1,M
BSAVE(I)=BETA(I)
BVAL(I)=BETA(I)
PIE(I)=PI(I)
LABEL(I)=IBASE(I)
IDEX(I)=I
INPROB(I)=I
SAVTOT(I)=TOT(I)
ZL(I)=BIG
JINM=JINM+ISTAT(I)
30 CONTINUE
DO 40 I=1,N
SIG(I)=SIGMA(I)
40 CONTINUE
C
LEVEL=0
DIRECT=.TRUE.
C
IF (JINM.GT.MININ) MININ=JINM
C
C INITIALIZE NUMBER OF FREE PARAMETERS IN MODEL
C
NFREE=M-JINM
C
ZL(M)=SAD
ZSAVE(M)=SAD
C
IF (MININ.EQ.M) RETURN
C
C NUMIN GIVES THE NUMBER OF PARAMETERS IN THE MODEL
C
NUMIN=M
C
C START THE MAIN LOOP
C
50 NUMIN=NUMIN-1
LEVEL=LEVEL+1
IF (NUMIN.LT.MININ.OR.NFREE.EQ.0) GO TO 120
K2=(M*(M+1)-(NUMIN+1)*(NUMIN+2))/2
KKK1=NUMIN+1
ZLOW=ZL(NUMIN)*POPT1
IF (DIRECT) GO TO 80
DO 60 I=1,KKK1
J=K2+I
INDEX(I)=I
IBASE(I)=LABEL(J)
INCOL(I)=INPROB(J)
PI(I)=PIE(J)
BETA(I)=BSAVE(J)
TOT(I)=SAVTOT(J)
60 CONTINUE
C
C LOAD IN THE BOUNDS FROM THE IMMEDIATE PREDECESSOR
C
K3=(M-KKK1)*N
DO 70 I=1,N
K3=K3+1
SIGMA(I)=SIG(K3)
70 CONTINUE
SAD=ZSAVE(KKK1)
C
C SET UP A NEW PROBLEM
C
80 CALL SETUP (X,KKK1,IFAULT,ISTAT,BETA,TOT,PI,REPP)
IF (IFAULT.GT.0) RETURN
NFREE=NFREE-1
C
C CALL PRIMAL L.P. CODE
C
CALL PHASE2 (X,Y,NUMIN,N,ITER,IFAULT,BETA,TOT,PI,REPP)
C
C SAVE THE SOLUTION DATA FOR LATER RECALL
C
K2=K2+KKK1+1
J=K2
DO 90 I=1,NUMIN
L=K2+I
PIE(L)=PI(I)
LABEL(L)=IBASE(I)
BSAVE(L)=BETA(I)
INPROB(L)=INCOL(I)
SAVTOT(L)=TOT(I)
90 CONTINUE
C
C SAVE THE OBJECTIVE VALUE
C
ZSAVE(NUMIN)=SAD
C
C SAVE THE NONBASIC BOUNDS
C
K3=(M-NUMIN)*N
DO 100 I=1,N
K3=K3+1
SIG(K3)=SIGMA(I)
100 CONTINUE
C
DIRECT=.TRUE.
C
C CHECK FOR A BETTER SOLUTION
C
IF (SADT.GE.ZLOW) GO TO 50
C
C SAVE THE BETTER SOLUTION
C
ZL(NUMIN)=SAD
K2=J
DO 110 I=1,NUMIN
L=K2+I
IDEX(L)=INCOL(I)
BVAL(L)=BETA(I)
110 CONTINUE
C
C GO TO TOP OF LOOP
C
GO TO 50
120 NUMIN=NUMIN+1
C
DIRECT=.FALSE.
C
C MUST WORK BACK UP THE TREE
C
LEVEL=LEVEL-1
130 IF (IPAR(LEVEL).GT.0) GO TO 140
IPAR(LEVEL)=-IPAR(LEVEL)
J=IPAR(LEVEL)
ISTAT(J)=1
NUMIN=NUMIN+1
GO TO 50
140 J=IPAR(LEVEL)
ISTAT(J)=0
NFREE=NFREE+1
LEVEL=LEVEL-1
IF (LEVEL.GT.0) GO TO 130
RETURN
C
END
SUBROUTINE SETUP (X,M,IFAULT,ISTAT,BETA,TOT,PI,REPP)
C
C THIS SUBROUTINE DETERMINES THE FORM OF THE CURRENT SUBPROBLEM
C
DIMENSION X(300,20), ISTAT(20), RHS(20), BETA(20), TOT(20), PI(20)
1, REPP(20)
C
COMMON /BLOCK/
- ACU,BIG,BSAVE(210),LU(20,20),PIE(210),SAD,SADT,SAVTOT(210),
-SIG(6000),SIGMA(300),ZLOW,ZSAVE(20),IBASE(20),IK,IK1,INCOL(20),
-INDEX(20),INPROB(210),IPAR(20),LABEL(210),LEVEL,DIRECT,INTL
C
INTEGER I,IBASE,IFAULT,IFIXI,II,IK,IK1,INCOL,INDEX,INPROB
INTEGER IOUT,IPAR,ISAVE,ISTAT,KKK,L,LABEL,LEAVE,LEVEL,LL,M
C
DOUBLE PRECISION ACU,BESIDE,BETA,BIG,BSAVE,LU,ONE,PENALT,PI,PIE
C++ REAL ACU,BESIDE,BETA,BIG,BSAVE,LU,ONE,PENALT,PI,PIE
DOUBLE PRECISION RATIO,REPP,RHO,RHS,SAD,SADT,SAVTOT,SIDE,SIG
C++ REAL RATIO,REPP,RHO,RHS,SAD,SADT,SAVTOT,SIDE,SIG
DOUBLE PRECISION SIGMA,SREPP,TEST,TOT,X,ZERO,ZLOW,ZSAVE
C++ REAL SIGMA,SREPP,TEST,TOT,X,ZERO,ZLOW,ZSAVE
C
LOGICAL DIRECT,INTL
C
DATA ONE/1.0D0/,ZERO/0.0D0/
C++ DATA ONE/1.0E0/,ZERO/0.0E0/
C
C IF THE IMMEDIATE PREDECESSOR IS IN MEMORY GO ON
C
IF (DIRECT) GO TO 10
C
C RECONSTRUCT THE LU DECOMPOSITION
C
KKK=1
CALL UPDATE (KKK,IFAULT,X,1,M)
IF (IFAULT.GT.0) RETURN
C
C RECALCULATE THE BASIC PI VALUES
C
KKK=0
C
CALL CALCPI (KKK,PI,TOT,M)
C
C DETERMINE PARAMETER TO LEAVE BASED ON PENALTY
C
10 PENALT=BIG
KKK=0
DO 20 L=1,M
RHS(L)=ZERO
20 CONTINUE
C
DO 40 I=1,M
II=INDEX(I)
LL=INCOL(II)
IF (ISTAT(LL).EQ.1) GO TO 40
C++ RHO=SIGN(1.,BETA(II))
RHO=DSIGN(1.D0,BETA(II))
RHS(II)=ONE
KKK=I-1
CALL CALCPI (KKK,REPP,RHS,M)
RHS(II)=ZERO
TEST=BIG
DO 30 L=1,M
C++ IF ( ABS(REPP(L)).LE.ACU) GO TO 30
IF (DABS(REPP(L)).LE.ACU) GO TO 30
C++ SREPP=SIGN(1.,REPP(L))
SREPP=DSIGN(1.D0,REPP(L))
C++ RATIO=ABS((1.-RHO*SREPP*PI(L))/REPP(L))
RATIO=DABS((1.-RHO*SREPP*PI(L))/REPP(L))
IF (RATIO.GE.TEST) GO TO 30
SIDE=RHO*SREPP
TEST=RATIO
IOUT=L
30 CONTINUE
C
C SEE IF THE PENALTY IS LESS
C
C++ IF (TEST*ABS(BETA(II)).GE.PENALT) GO TO 40
IF (TEST*DABS(BETA(II)).GE.PENALT) GO TO 40
LEAVE=IOUT
IFIXI=II
BESIDE=SIDE
C++ PENALT=TEST*ABS(BETA(II))
PENALT=TEST*DABS(BETA(II))
40 CONTINUE
C
C UPDATE THE OBJECTIVE VALUE
C
SAD=SAD+PENALT
C
C PARAMETER INCOL(IFIXI) WILL LEAVE THE MODEL
C PI(IBASE(LEAVE)) WILL LEAVE THE BASIS AND GO TO BOUND
C
C SWITCH UNWANTED PARAMETER AND PI TO THE END OF LIST
C
ISAVE=INCOL(IFIXI)
INCOL(IFIXI)=INCOL(M)
INCOL(M)=ISAVE
TOT(IFIXI)=TOT(M)
C
C LABEL THIS NODE AND FLAG THE OUTGOING PARAMETER
C
IPAR(LEVEL)=-ISAVE
ISTAT(ISAVE)=-1
C
ISAVE=IBASE(LEAVE)
IBASE(LEAVE)=IBASE(M)
IBASE(M)=ISAVE
C
C PLACE THE OUTGOING PI AT THE CORRECT BOUND
C
SIGMA(ISAVE)=BESIDE
C
C FORM A NEW BASIS
C
M=M-1
C
DO 50 I=1,M
50 INDEX(I)=I
C
KKK=1
CALL UPDATE (KKK,IFAULT,X,1,M)
IF (IFAULT.GT.0) RETURN
C
C UPDATE THE TOT ARRAY
C
DO 60 I=1,M
LL=INCOL(I)
TOT(I)=TOT(I)-BESIDE*X(ISAVE,LL)
60 CONTINUE
C
RETURN
C
END
SUBROUTINE PHASE2 (X,Y,M,N,ITER,IFAULT,BETA,TOT,PI,REPP)
C
C SOLVES LP WITH PRIMAL ALGORITHM
C
DIMENSION X(300,20), Y(300), BETA(20), TOT(20), PI(20), REPP(20),
1RHS(20),Z(10),KAN(10)
C
COMMON /BLOCK/
- ACU,BIG,BSAVE(210),LU(20,20),PIE(210),SAD,SADT,SAVTOT(210),
-SIG(6000),SIGMA(300),ZLOW,ZSAVE(20),IBASE(20),IK,IK1,INCOL(20),
-INDEX(20),INPROB(210),IPAR(20),LABEL(210),LEVEL,DIRECT,INTL
C
INTEGER I,IBASE,IFAULT,IIN,IK,IK1,INCOL,INDEX,INPROB,IOUT,IPAR
INTEGER ITER,J,K,KAN,KKK,KROW,L,LABEL,LEVEL,LL,M,N
C
DOUBLE PRECISION ACU,BEST,BETA,BIG,BSAVE,LU,ONE,PI,PIE,RATIO,REPP
C++ REAL ACU,BEST,BETA,BIG,BSAVE,LU,ONE,PI,PIE,RATIO,REPP
DOUBLE PRECISION RHO,RHS,SAD,SADT,SAVTOT,SIG,SIGMA,SREPP,TEST
C++ REAL RHO,RHS,SAD,SADT,SAVTOT,SIG,SIGMA,SREPP,TEST
DOUBLE PRECISION TOT,X,XTWO,Y,Z,ZERO,ZC,ZLOW,ZSAVE,ZZ
C++ REAL TOT,X,XTWO,Y,Z,ZERO,ZC,ZLOW,ZSAVE,ZZ
C
LOGICAL DIRECT,INTL
C
DATA ZERO/0.0D0/,ONE/1.0D0/
C++ DATA ZERO/0.0E0/,ONE/1.0E0/
C
C Z(I) STORES THE REDUCED COST VALUES ON THE LIST
C KAN(I) STORES THE ROW NUMBER ON THE LIST
C IK IS THE LENGTH OF THE LIST
C
C CALCULATE THE SIMPLEX MULTIPLIERS
C
DO 10 I=1,M
K=IBASE(I)
10 PI(I)=Y(K)
KKK=0
CALL CALBET (KKK,BETA,PI,M)
C
C RECALCULATE THE BASIC PI VALUES
C
KKK=0
CALL CALCPI (KKK,PI,TOT,M)
C
C STORE THE OBJECTIVE VALUE USED FOR TERMINATION CHECK
C
SADT=SAD
IF (SAD.GE.ZLOW) RETURN
C
C GENERATE THE CANDIDATE LIST
C
20 DO 30 I=1,IK
Z(I)=ACU
30 CONTINUE
C
C CALCULATE THE REDUCED COSTS
C
DO 70 J=1,N
IF (SIGMA(J).EQ.ZERO) GO TO 70
ZC=ZERO
DO 40 I=1,M
LL=INCOL(I)
40 ZC=ZC+BETA(I)*X(J,LL)
ZZ=(ZC-Y(J))*SIGMA(J)
IF (ZZ.LE.Z(IK)) GO TO 70
C
C RANK THE VALUES IN Z(I), IN DESCENDING ORDER
C
DO 60 K=1,IK1
IF (ZZ.LE.Z(K)) GO TO 60
I=IK
50 Z(I)=Z(I-1)
KAN(I)=KAN(I-1)
I=I-1
IF (I.GT.K) GO TO 50
Z(K)=ZZ
KAN(K)=J
GO TO 70
60 CONTINUE
Z(IK)=ZZ
KAN(IK)=J
70 CONTINUE
IF (Z(1).LE.ACU) RETURN
KROW=1
ZZ=Z(1)
IIN=KAN(1)
ZC=ZZ*SIGMA(IIN)
GO TO 100
80 KROW=KROW+1
IF (KROW.GT.IK) GO TO 20
IF (Z(KROW).LE.ACU) GO TO 20
C
C DETERMINE THE ENTERING VARIABLE FROM THE CANDIDATE LIST
C
IIN=KAN(KROW)
ZC=ZERO
DO 90 I=1,M
LL=INCOL(I)
90 ZC=ZC+BETA(I)*X(IIN,LL)
ZC=ZC-Y(IIN)
ZZ=ZC*SIGMA(IIN)
IF (ZZ.LE.ACU) GO TO 80
100 BEST=ZZ
RHO=SIGMA(IIN)
C
C RHO: SIGN OF THE INCOMING REDUCED COST
C IIN: INCOMING CANDIDATE
C
C FIND THE REPRESENTATION OF THE ENTERING VARIABLE
C
DO 110 I=1,M
LL=INCOL(I)
RHS(I)=X(IIN,LL)
110 CONTINUE
KKK=0
CALL CALCPI (KKK,REPP,RHS,M)
C
C CALCULATE THE MIN RATIO TEST TO DETERMINE THE LEAVING VARIABLE
C
TEST=2.0
DO 120 I=1,M
C++ IF (ABS(REPP(I)).LE.ACU) GO TO 120
IF (DABS(REPP(I)).LE.ACU) GO TO 120
C++ SREPP=SIGN(1.0,REPP(I))
SREPP=DSIGN(1.0D0,REPP(I))
C++ RATIO=ABS((1.-RHO*SREPP*PI(I))/REPP(I))
RATIO=DABS((1.-RHO*SREPP*PI(I))/REPP(I))
IF (RATIO.GE.TEST) GO TO 120
TEST=RATIO
IOUT=I
120 CONTINUE
C
C PERFORM FATHOMING TEST BEFORE PIVOT
C
SADT=SAD+TEST*ZZ
IF (SADT.GE.ZLOW) RETURN
SAD=SADT
C
DO 130 I=1,M
130 PI(I)=PI(I)+TEST*REPP(I)*RHO
IF (TEST.LT.2.0) GO TO 150
C
SIGMA(IIN)=-SIGMA(IIN)
XTWO=SIGMA(IIN)+SIGMA(IIN)
DO 140 I=1,M
LL=INCOL(I)
TOT(I)=TOT(I)-XTWO*X(IIN,LL)
140 CONTINUE
GO TO 80
150 K=IBASE(IOUT)
C++ SIGMA(K)=SIGN(1.0,PI(IOUT))
SIGMA(K)=DSIGN(1.0D0,PI(IOUT))
PI(IOUT)=SIGMA(IIN)*(1.0-TEST)
IBASE(IOUT)=IIN
C
DO 160 I=1,M
LL=INCOL(I)
TOT(I)=TOT(I)+SIGMA(IIN)*X(IIN,LL)-SIGMA(K)*X(K,LL)
160 CONTINUE
C
KKK=IOUT
CALL UPDATE (IOUT,IFAULT,X,1,M)
C
DO 170 I=1,M
RHS(I)=ZERO
170 CONTINUE
RHS(KKK)=ONE
KKK=KKK-1
CALL CALBET (KKK,REPP,RHS,M)
DO 180 L=1,M
180 BETA(L)=BETA(L)-ZC*REPP(L)
SIGMA(IIN)=ZERO
ITER=ITER+1
GO TO 80
C
END
SUBROUTINE L1NORM (X,Y,N,M,IFAULT,ITER,BETA,PI,TOT)
C
C THIS ROUTINE SOLVES THE INITIAL REGRESSION PROBLEM WITH ALL
C THE PARAMETERS INCLUDED IN THE MODEL
C
DIMENSION X(300,20), Y(300), D(300), DELTA(300), INEXT(300), PRICE
1(300), RHS(20), BETA(20), TOT(20), PI(20)
C
C
COMMON /BLOCK/
- ACU,BIG,BSAVE(210),LU(20,20),PIE(210),SAD,SADT,SAVTOT(210),
-SIG(6000),SIGMA(300),ZLOW,ZSAVE(20),IBASE(20),IK,IK1,INCOL(20),
-INDEX(20),INPROB(210),IPAR(20),LABEL(210),LEVEL,DIRECT,INTL
C
INTEGER I,IBASE,IFAULT,IIN,IK,IK1,INCOL,INEXT,INPROB,IOUT,IPAR
INTEGER INDEX,IPT,ITER,J,K,K1,KKK,KOUNT,L,LABEL,LEVEL,M,M1,N
C
DOUBLE PRECISION ACU,AHALF,AONE,BETA,BIG,BONE,BSAVE,D,DELTA,LU
C++ REAL ACU,AHALF,AONE,BETA,BIG,BONE,BSAVE,D,DELTA,LU
DOUBLE PRECISION ONE,PI,PIE,PRICE,RATIO,RHO,RHS,SAD,SADT,SAVTOT
C++ REAL ONE,PI,PIE,PRICE,RATIO,RHO,RHS,SAD,SADT,SAVTOT
DOUBLE PRECISION SIG,SIGMA,SUBT,SUM,T,TEST,TOT,VAL,X,Y,ZERO,ZLOW
C++ REAL SIG,SIGMA,SUBT,SUM,T,TEST,TOT,VAL,X,Y,ZERO,ZLOW
DOUBLE PRECISION ZSAVE,ZZZ
C++ REAL ZSAVE,ZZZ
C
LOGICAL DIRECT,INTL
C
DATA ONE/1.0D0/,ZERO/0.0D0/
C++ DATA ONE/1.0E0/,ZERO/0.0E0/
C
C INITIAL SETTINGS
C
ITER=0
AONE=ONE+ACU
AHALF=.5+ACU
BONE=ONE-ACU
M1=M-1
C
C SET UP INITIAL LU DECOMPOSITION
C
INTL=.TRUE.
K=1
CALL UPDATE (K,IFAULT,X,N,M)
IF (IFAULT.NE.0) RETURN
INTL=.FALSE.
DO 10 I=1,M
K1=IBASE(I)
RHS(I)=Y(K1)
10 CONTINUE
KKK=0
CALL CALBET (KKK,BETA,RHS,M)
C
C CALCULATE INITIAL D, TOT AND SIGMA VECTORS
C
DO 30 J=1,N
VAL=ZERO
DO 20 I=1,M
VAL=VAL+BETA(I)*X(J,I)
20 CONTINUE
D(J)=Y(J)-VAL
C++ SIGMA(J)=SIGN(1.,D(J))
SIGMA(J)=DSIGN(1.D0,D(J))
30 CONTINUE
DO 40 J=1,M
RHS(J)=ZERO
KKK=IBASE(J)
SIGMA(KKK)=ZERO
40 CONTINUE
DO 60 J=1,N
DO 50 I=1,M
50 TOT(I)=TOT(I)-SIGMA(J)*X(J,I)
60 CONTINUE
C
C MAIN ITERATIVE LOOP BEGINS
C
70 CONTINUE
C
KKK=0
CALL CALCPI (KKK,PI,TOT,M)
C
T=AONE
K=0
DO 80 J=1,M
C++ IF (ABS(PI(J)).LT.T) GO TO 80
IF (DABS(PI(J)).LT.T) GO TO 80
K=J
C++ T=ABS(PI(J))
T=DABS(PI(J))
C++ RHO=-SIGN(1.,PI(J))
RHO=-DSIGN(1.D0,PI(J))
80 CONTINUE
IF (K.EQ.0) GO TO 220
KKK=K-1
RHS(K)=ONE
CALL CALBET (KKK,BETA,RHS,M)
RHS(K)=ZERO
DO 100 I=1,N
DELTA(I)=ZERO
IF (SIGMA(I).EQ.ZERO) GO TO 100
DO 90 J=1,M
DELTA(I)=DELTA(I)+BETA(J)*X(I,J)
90 CONTINUE
DELTA(I)=RHO*DELTA(I)
100 CONTINUE
C
C PERFORM PARTIAL SORT OF RATIOS
C
T=T*.5
KOUNT=0
RATIO=BIG
SUM=AHALF
SUBT=ZERO
DO 160 I=1,N
IF (DELTA(I)*SIGMA(I).LE.ACU) GO TO 160
TEST=D(I)/DELTA(I)
IF (TEST.GE.RATIO) GO TO 160
C++ SUM=SUM+ABS(DELTA(I))
SUM=SUM+DABS(DELTA(I))
IF (SUM-SUBT.GE.T) GO TO 110
C
C INSERT I IN LIST
C
KOUNT=KOUNT+1
PRICE(KOUNT)=TEST
INEXT(KOUNT)=I
GO TO 160
C
C UPDATE SUM AND KICK IIN OUT OF THE LIST
C
110 SUM=SUM-SUBT
RATIO=TEST
IPT=0
KKK=0
C
C IDENTIFY A NEW IIN
C
120 KKK=KKK+1
IF (KKK.GT.KOUNT) GO TO 130
IF (PRICE(KKK).LE.RATIO) GO TO 120
RATIO=PRICE(KKK)
IPT=KKK
GO TO 120
130 IF (IPT.EQ.0) GO TO 150
C
C SWITCH VALUES
C
KKK=INEXT(IPT)
C++ SUBT=ABS(DELTA(KKK))
SUBT=DABS(DELTA(KKK))
IF (SUM-SUBT.LT.T) GO TO 140
PRICE(IPT)=PRICE(KOUNT)
INEXT(IPT)=INEXT(KOUNT)
KOUNT=KOUNT-1
GO TO 110
140 IIN=INEXT(IPT)
INEXT(IPT)=I
PRICE(IPT)=TEST
GO TO 160
150 IIN=I
C++ SUBT=ABS(DELTA(I))
SUBT=DABS(DELTA(I))
160 CONTINUE
C
C UPDATE BASIC INDICATORS
C
IF (KOUNT.EQ.0) GO TO 190
DO 180 J=1,KOUNT
KKK=INEXT(J)
ZZZ=SIGMA(KKK)
DO 170 L=1,M
SUBT=ZZZ*X(KKK,L)
TOT(L)=TOT(L)+SUBT+SUBT
170 CONTINUE
SIGMA(KKK)=-SIGMA(KKK)
180 CONTINUE
190 CONTINUE
IOUT=IBASE(K)
DELTA(IOUT)=RHO
DO 200 L=1,M
200 TOT(L)=TOT(L)+RHO*X(IOUT,L)+SIGMA(IIN)*X(IIN,L)
SIGMA(IOUT)=-RHO
IBASE(K)=IIN
CALL UPDATE (K,IFAULT,X,1,M)
DO 210 J=1,N
210 D(J)=D(J)-RATIO*DELTA(J)
SIGMA(IIN)=ZERO
ITER=ITER+1
GO TO 70
C
C CALCULATE OPTIMAL BETA AND SUM OF ABSOLUTE DEVIATIONS
C
220 DO 230 I=1,M
K1=IBASE(I)
RHS(I)=Y(K1)
230 CONTINUE
KKK=0
CALL CALBET (KKK,BETA,RHS,M)
SAD=ZERO
DO 240 J=1,N
C++ 240 SAD=SAD+ABS(D(J))
240 SAD=SAD+DABS(D(J))
RETURN
C
END
SUBROUTINE CALBET (KKK,BETA,RHS,M)
C
C SUBROUTINE CALBET FOR BACK-SOLVING SYSTEM OF EQUATIONS
C
DIMENSION BETA(20), RHS(20)
C
COMMON /BLOCK/
- ACU,BIG,BSAVE(210),LU(20,20),PIE(210),SAD,SADT,SAVTOT(210),
-SIG(6000),SIGMA(300),ZLOW,ZSAVE(20),IBASE(20),IK,IK1,INCOL(20),
-INDEX(20),INPROB(210),IPAR(20),LABEL(210),LEVEL,DIRECT,INTL
C
INTEGER I,IBASE,II,II1,IK,IK1,INCOL,INDEX,INPROB,IPAR,K,K1,K2
INTEGER KK,KKK,KKK1,LABEL,LEVEL,M,M1
C
DOUBLE PRECISION ACU,BETA,BIG,BSAVE,LU,PIE,RHS,SAD,SADT,SAVTOT
C++ REAL ACU,BETA,BIG,BSAVE,LU,PIE,RHS,SAD,SADT,SAVTOT
DOUBLE PRECISION SIG,SIGMA,ZERO,ZLOW,ZSAVE
C++ REAL SIG,SIGMA,ZERO,ZLOW,ZSAVE
C
LOGICAL DIRECT,INTL
C
DATA ZERO/0.0D0/
C++ DATA ZERO/0.0E0/
C
IF (M.GT.1) GO TO 10
K=INDEX(1)
BETA(K)=RHS(1)/LU(K,1)
RETURN
10 CONTINUE
M1=M-1
IF (KKK.EQ.0) GO TO 30
DO 20 I=1,KKK
K=INDEX(I)
BETA(K)=ZERO
20 CONTINUE
30 KKK=KKK+1
K=INDEX(KKK)
BETA(K)=RHS(KKK)/LU(K,KKK)
IF (KKK.EQ.M) GO TO 60
KKK1=KKK+1
DO 50 II=KKK1,M
K=INDEX(II)
BETA(K)=RHS(II)
II1=II-1
DO 40 I=KKK,II1
KK=INDEX(I)
BETA(K)=BETA(K)-LU(KK,II)*BETA(KK)
40 CONTINUE
BETA(K)=BETA(K)/LU(K,II)
50 CONTINUE
60 CONTINUE
DO 70 II=1,M1
K1=M-II
K=INDEX(K1)
DO 70 I=1,II
KK=M-I+1
K2=INDEX(KK)
BETA(K)=BETA(K)-LU(K2,K1)*BETA(K2)
70 CONTINUE
RETURN
C
END
SUBROUTINE UPDATE (KKK,IFAULT,X,N,M)
C
C SUBROUTINE UPDATE FOR UPDATING LU DECOMPOSITION MATRIX
C
COMMON /BLOCK/
- ACU,BIG,BSAVE(210),LU(20,20),PIE(210),SAD,SADT,SAVTOT(210),
-SIG(6000),SIGMA(300),ZLOW,ZSAVE(20),IBASE(20),IK,IK1,INCOL(20),
-INDEX(20),INPROB(210),IPAR(20),LABEL(210),LEVEL,DIRECT,INTL
C
INTEGER I,IBASE,ICOL,IFAULT,II,II1,IK,IK1,INCOL,INDEX,INPROB
INTEGER IPAR,IROW,ISAVE,J,K,KK,KKK,LABEL,LEVEL,LL,M,N
C
DOUBLE PRECISION ACU,BIG,BSAVE,LU,PIE,PIVOT,SAD,SADT,SAVTOT
C++ REAL ACU,BIG,BSAVE,LU,PIE,PIVOT,SAD,SADT,SAVTOT
DOUBLE PRECISION SIG,SIGMA,SUBT,X,ZLOW,ZSAVE
C++ REAL SIG,SIGMA,SUBT,X,ZLOW,ZSAVE
C
LOGICAL DIRECT,INTL
C
DIMENSION X(300,20)
C
IROW=0
DO 90 II=KKK,M
IF (INTL) GO TO 10
IROW=IBASE(II)
GO TO 20
10 IROW=IROW+1
IBASE(II)=IROW
IF (IROW.GT.N) IFAULT=1
IF (IFAULT.NE.0) RETURN
20 DO 30 I=1,M
LL=INCOL(I)
LU(I,II)=X(IROW,LL)
30 CONTINUE
C
C SET UP REPRESENTATION OF INCOMING ROW
C
IF (II.EQ.1) GO TO 60
II1=II-1
DO 50 ICOL=1,II1
K=INDEX(ICOL)
SUBT=LU(K,II)
J=ICOL+1
DO 40 I=J,M
K=INDEX(I)
LU(K,II)=LU(K,II)-SUBT*LU(K,ICOL)
40 CONTINUE
50 CONTINUE
C
C FIND MAXIMUM ENTRY
C
60 PIVOT=ACU
KK=0
DO 70 I=II,M
K=INDEX(I)
IF (DABS(LU(K,II)).LE.PIVOT) GO TO 70
C++ IF (ABS(LU(K,II)).LE.PIVOT) GO TO 70
C++ PIVOT=ABS(LU(K,II))
PIVOT=DABS(LU(K,II))
KK=I
70 CONTINUE
IF (KK.EQ.0) GO TO 10
C
C SWITCH ORDER
C
ISAVE=INDEX(KK)
INDEX(KK)=INDEX(II)
INDEX(II)=ISAVE
C
C PUT IN COLUMNS OF LU ONE AT A TIME
C
IF (II.EQ.M) GO TO 90
J=II+1
DO 80 I=J,M
K=INDEX(I)
LU(K,II)=LU(K,II)/LU(ISAVE,II)
80 CONTINUE
90 CONTINUE
KKK=IROW
RETURN
C
END
SUBROUTINE CALCPI (KKK,PI,TOT,M)
C
C THIS ROUTINE FORWARD SOLVES A LINEAR SYSTEM
C
DIMENSION TOT(20), PI(20)
C
COMMON /BLOCK/
- ACU,BIG,BSAVE(210),LU(20,20),PIE(210),SAD,SADT,SAVTOT(210),
-SIG(6000),SIGMA(300),ZLOW,ZSAVE(20),IBASE(20),IK,IK1,INCOL(20),
-INDEX(20),INPROB(210),IPAR(20),LABEL(210),LEVEL,DIRECT,INTL
C
INTEGER I,IBASE,II,II1,IK,IK1,INCOL,INDEX,INPROB,IPAR,K,K1,K2,KKK
INTEGER KKK1,LABEL,LEVEL,M,M1
C
DOUBLE PRECISION ACU,BIG,BSAVE,LU,PI,PIE,SAD,SADT,SAVTOT,SIG
C++ REAL ACU,BIG,BSAVE,LU,PI,PIE,SAD,SADT,SAVTOT,SIG
DOUBLE PRECISION SIGMA,TOT,ZERO,ZLOW,ZSAVE
C++ REAL SIGMA,TOT,ZERO,ZLOW,ZSAVE
C
LOGICAL DIRECT,INTL
C
C++ DATA ZERO/0.0E0/
DATA ZERO/0.0D0/
M1=M-1
C
IF (M.GT.1) GO TO 10
K=INDEX(M)
PI(M)=TOT(K)/LU(K,M)
RETURN
10 IF (KKK.EQ.0) GO TO 30
DO 20 I=1,KKK
PI(I)=ZERO
20 CONTINUE
30 CONTINUE
KKK=KKK+1
K=INDEX(KKK)
PI(KKK)=TOT(K)
IF (KKK.EQ.M) GO TO 60
KKK1=KKK+1
DO 50 II=KKK1,M
K=INDEX(II)
PI(II)=TOT(K)
II1=II-1
DO 40 I=KKK,II1
40 PI(II)=PI(II)-LU(K,I)*PI(I)
50 CONTINUE
60 CONTINUE
K=INDEX(M)
PI(M)=PI(M)/LU(K,M)
DO 80 II=1,M1
K1=M-II
K=INDEX(K1)
DO 70 I=1,II
K2=M-I+1
PI(K1)=PI(K1)-LU(K,K2)*PI(K2)
70 CONTINUE
PI(K1)=PI(K1)/LU(K,K1)
80 CONTINUE
RETURN
END
DOUBLE PRECISION FUNCTION D1MACH(I)
C***BEGIN PROLOGUE D1MACH
C***DATE WRITTEN 750101 (YYMMDD)
C***REVISION DATE 820801 (YYMMDD)
C***CATEGORY NO. R1
C***KEYWORDS MACHINE CONSTANTS
C***AUTHOR FOX, P. A., (BELL LABS)
C HALL, A. D., (BELL LABS)
C SCHRYER, N. L., (BELL LABS)
C***PURPOSE RETURNS DOUBLE PRECISION MACHINE DEPENDENT CONSTANTS
C***DESCRIPTION
C D1MACH CAN BE USED TO OBTAIN MACHINE-DEPENDENT PARAMETERS
C FOR THE LOCAL MACHINE ENVIRONMENT. IT IS A FUNCTION
C SUBPROGRAM WITH ONE (INPUT) ARGUMENT, AND CAN BE CALLED
C AS FOLLOWS, FOR EXAMPLE
C
C D = D1MACH(I)
C
C WHERE I=1,...,5. THE (OUTPUT) VALUE OF D ABOVE IS
C DETERMINED BY THE (INPUT) VALUE OF I. THE RESULTS FOR
C VARIOUS VALUES OF I ARE DISCUSSED BELOW.
C
C DOUBLE-PRECISION MACHINE CONSTANTS
C D1MACH( 1) = B**(EMIN-1), THE SMALLEST POSITIVE MAGNITUDE.
C D1MACH( 2) = B**EMAX*(1 - B**(-T)), THE LARGEST MAGNITUDE.
C D1MACH( 3) = B**(-T), THE SMALLEST RELATIVE SPACING.
C D1MACH( 4) = B**(1-T), THE LARGEST RELATIVE SPACING.
C D1MACH( 5) = LOG10(B)
C***REFERENCES FOX P.A., HALL A.D., SCHRYER N.L.,*FRAMEWORK FOR A
C PORTABLE LIBRARY*, ACM TRANSACTIONS ON MATHEMATICAL
C SOFTWARE, VOL. 4, NO. 2, JUNE 1978, PP. 177-188.
C***ROUTINES CALLED XERROR
C***END PROLOGUE D1MACH
C
INTEGER SMALL(4)
INTEGER LARGE(4)
INTEGER RIGHT(4)
INTEGER DIVER(4)
INTEGER LOG10(4)
C
DOUBLE PRECISION DMACH(5)
C
EQUIVALENCE (DMACH(1),SMALL(1))
EQUIVALENCE (DMACH(2),LARGE(1))
EQUIVALENCE (DMACH(3),RIGHT(1))
EQUIVALENCE (DMACH(4),DIVER(1))
EQUIVALENCE (DMACH(5),LOG10(1))
C
C MACHINE CONSTANTS FOR THE BURROUGHS 1700 SYSTEM.
C
C DATA SMALL(1) / ZC00800000 /
C DATA SMALL(2) / Z000000000 /
C
C DATA LARGE(1) / ZDFFFFFFFF /
C DATA LARGE(2) / ZFFFFFFFFF /
C
C DATA RIGHT(1) / ZCC5800000 /
C DATA RIGHT(2) / Z000000000 /
C
C DATA DIVER(1) / ZCC6800000 /
C DATA DIVER(2) / Z000000000 /
C
C DATA LOG10(1) / ZD00E730E7 /
C DATA LOG10(2) / ZC77800DC0 /
C
C MACHINE CONSTANTS FOR THE BURROUGHS 5700 SYSTEM.
C
C DATA SMALL(1) / O1771000000000000 /
C DATA SMALL(2) / O0000000000000000 /
C
C DATA LARGE(1) / O0777777777777777 /
C DATA LARGE(2) / O0007777777777777 /
C
C DATA RIGHT(1) / O1461000000000000 /
C DATA RIGHT(2) / O0000000000000000 /
C
C DATA DIVER(1) / O1451000000000000 /
C DATA DIVER(2) / O0000000000000000 /
C
C DATA LOG10(1) / O1157163034761674 /
C DATA LOG10(2) / O0006677466732724 /
C
C MACHINE CONSTANTS FOR THE BURROUGHS 6700/7700 SYSTEMS.
C
C DATA SMALL(1) / O1771000000000000 /
C DATA SMALL(2) / O7770000000000000 /
C
C DATA LARGE(1) / O0777777777777777 /
C DATA LARGE(2) / O7777777777777777 /
C
C DATA RIGHT(1) / O1461000000000000 /
C DATA RIGHT(2) / O0000000000000000 /
C
C DATA DIVER(1) / O1451000000000000 /
C DATA DIVER(2) / O0000000000000000 /
C
C DATA LOG10(1) / O1157163034761674 /
C DATA LOG10(2) / O0006677466732724 /
C
C MACHINE CONSTANTS FOR THE CDC 6000/7000 SERIES.
C
C DATA SMALL(1) / 00604000000000000000B /
C DATA SMALL(2) / 00000000000000000000B /
C
C DATA LARGE(1) / 37767777777777777777B /
C DATA LARGE(2) / 37167777777777777777B /
C
C DATA RIGHT(1) / 15604000000000000000B /
C DATA RIGHT(2) / 15000000000000000000B /
C
C DATA DIVER(1) / 15614000000000000000B /
C DATA DIVER(2) / 15010000000000000000B /
C
C DATA LOG10(1) / 17164642023241175717B /
C DATA LOG10(2) / 16367571421742254654B /
C
C MACHINE CONSTANTS FOR THE CRAY 1
C
C DATA SMALL(1) / 200004000000000000000B /
C DATA SMALL(2) / 00000000000000000000B /
C
C DATA LARGE(1) / 577777777777777777777B /
C DATA LARGE(2) / 000007777777777777777B /
C
C DATA RIGHT(1) / 377214000000000000000B /
C DATA RIGHT(2) / 000000000000000000000B /
C
C DATA DIVER(1) / 377224000000000000000B /
C DATA DIVER(2) / 000000000000000000000B /
C
C DATA LOG10(1) / 377774642023241175717B /
C DATA LOG10(2) / 000007571421742254654B /
C
C MACHINE CONSTANTS FOR THE DATA GENERAL ECLIPSE S/200
C
C NOTE - IT MAY BE APPROPRIATE TO INCLUDE THE FOLLOWING CARD -
C STATIC DMACH(5)
C
C DATA SMALL/20K,3*0/,LARGE/77777K,3*177777K/
C DATA RIGHT/31420K,3*0/,DIVER/32020K,3*0/
C DATA LOG10/40423K,42023K,50237K,74776K/
C
C MACHINE CONSTANTS FOR THE HARRIS 220
C
C DATA SMALL(1),SMALL(2) / [20000000, [00000201 /
C DATA LARGE(1),LARGE(2) / [37777777, [37777577 /
C DATA RIGHT(1),RIGHT(2) / [20000000, [00000333 /
C DATA DIVER(1),DIVER(2) / [20000000, [00000334 /
C DATA LOG10(1),LOG10(2) / [23210115, [10237777 /
C
C MACHINE CONSTANTS FOR THE HONEYWELL 600/6000 SERIES.
C
C DATA SMALL(1),SMALL(2) / O402400000000, O000000000000 /
C DATA LARGE(1),LARGE(2) / O376777777777, O777777777777 /
C DATA RIGHT(1),RIGHT(2) / O604400000000, O000000000000 /
C DATA DIVER(1),DIVER(2) / O606400000000, O000000000000 /
C DATA LOG10(1),LOG10(2) / O776464202324, O117571775714 /
C
C MACHINE CONSTANTS FOR THE HP 2100
C THREE WORD DOUBLE PRECISION OPTION WITH FTN4
C
C DATA SMALL(1), SMALL(2), SMALL(3) / 40000B, 0, 1 /
C DATA LARGE(1), LARGE(2), LARGE(3) / 77777B, 177777B, 177776B /
C DATA RIGHT(1), RIGHT(2), RIGHT(3) / 40000B, 0, 265B /
C DATA DIVER(1), DIVER(2), DIVER(3) / 40000B, 0, 276B /
C DATA LOG10(1), LOG10(2), LOG10(3) / 46420B, 46502B, 77777B /
C
C
C MACHINE CONSTANTS FOR THE HP 2100
C FOUR WORD DOUBLE PRECISION OPTION WITH FTN4
C
C DATA SMALL(1), SMALL(2) / 40000B, 0 /
C DATA SMALL(3), SMALL(4) / 0, 1 /
C DATA LARGE(1), LARGE(2) / 77777B, 177777B /
C DATA LARGE(3), LARGE(4) / 177777B, 177776B /
C DATA RIGHT(1), RIGHT(2) / 40000B, 0 /
C DATA RIGHT(3), RIGHT(4) / 0, 225B /
C DATA DIVER(1), DIVER(2) / 40000B, 0 /
C DATA DIVER(3), DIVER(4) / 0, 227B /
C DATA LOG10(1), LOG10(2) / 46420B, 46502B /
C DATA LOG10(3), LOG10(4) / 76747B, 176377B /
C
C
C MACHINE CONSTANTS FOR THE IBM 360/370 SERIES,
C THE XEROX SIGMA 5/7/9, THE SEL SYSTEMS 85/86, AND
C THE PERKIN ELMER (INTERDATA) 7/32.
C
C DATA SMALL(1),SMALL(2) / Z00100000, Z00000000 /
C DATA LARGE(1),LARGE(2) / Z7FFFFFFF, ZFFFFFFFF /
C DATA RIGHT(1),RIGHT(2) / Z33100000, Z00000000 /
C DATA DIVER(1),DIVER(2) / Z34100000, Z00000000 /
C DATA LOG10(1),LOG10(2) / Z41134413, Z509F79FF /
C
C MACHINE CONSTANTS FOR THE PDP-10 (KA PROCESSOR).
C
C DATA SMALL(1),SMALL(2) / "033400000000, "000000000000 /
C DATA LARGE(1),LARGE(2) / "377777777777, "344777777777 /
C DATA RIGHT(1),RIGHT(2) / "113400000000, "000000000000 /
C DATA DIVER(1),DIVER(2) / "114400000000, "000000000000 /
C DATA LOG10(1),LOG10(2) / "177464202324, "144117571776 /
C
C MACHINE CONSTANTS FOR THE PDP-10 (KI PROCESSOR).
C
C DATA SMALL(1),SMALL(2) / "000400000000, "000000000000 /
C DATA LARGE(1),LARGE(2) / "377777777777, "377777777777 /
C DATA RIGHT(1),RIGHT(2) / "103400000000, "000000000000 /
C DATA DIVER(1),DIVER(2) / "104400000000, "000000000000 /
C DATA LOG10(1),LOG10(2) / "177464202324, "476747767461 /
C
C MACHINE CONSTANTS FOR PDP-11 FORTRAN SUPPORTING
C 32-BIT INTEGERS (EXPRESSED IN INTEGER AND OCTAL).
C
C DATA SMALL(1),SMALL(2) / 8388608, 0 /
C DATA LARGE(1),LARGE(2) / 2147483647, -1 /
C DATA RIGHT(1),RIGHT(2) / 612368384, 0 /
C DATA DIVER(1),DIVER(2) / 620756992, 0 /
C DATA LOG10(1),LOG10(2) / 1067065498, -2063872008 /
C
C DATA SMALL(1),SMALL(2) / O00040000000, O00000000000 /
C DATA LARGE(1),LARGE(2) / O17777777777, O37777777777 /
C DATA RIGHT(1),RIGHT(2) / O04440000000, O00000000000 /
C DATA DIVER(1),DIVER(2) / O04500000000, O00000000000 /
C DATA LOG10(1),LOG10(2) / O07746420232, O20476747770 /
C
C MACHINE CONSTANTS FOR PDP-11 FORTRAN SUPPORTING
C 16-BIT INTEGERS (EXPRESSED IN INTEGER AND OCTAL).
C
C DATA SMALL(1),SMALL(2) / 128, 0 /
C DATA SMALL(3),SMALL(4) / 0, 0 /
C
C DATA LARGE(1),LARGE(2) / 32767, -1 /
C DATA LARGE(3),LARGE(4) / -1, -1 /
C
C DATA RIGHT(1),RIGHT(2) / 9344, 0 /
C DATA RIGHT(3),RIGHT(4) / 0, 0 /
C
C DATA DIVER(1),DIVER(2) / 9472, 0 /
C DATA DIVER(3),DIVER(4) / 0, 0 /
C
C DATA LOG10(1),LOG10(2) / 16282, 8346 /
C DATA LOG10(3),LOG10(4) / -31493, -12296 /
C
C DATA SMALL(1),SMALL(2) / O000200, O000000 /
C DATA SMALL(3),SMALL(4) / O000000, O000000 /
C
C DATA LARGE(1),LARGE(2) / O077777, O177777 /
C DATA LARGE(3),LARGE(4) / O177777, O177777 /
C
C DATA RIGHT(1),RIGHT(2) / O022200, O000000 /
C DATA RIGHT(3),RIGHT(4) / O000000, O000000 /
C
C DATA DIVER(1),DIVER(2) / O022400, O000000 /
C DATA DIVER(3),DIVER(4) / O000000, O000000 /
C
C DATA LOG10(1),LOG10(2) / O037632, O020232 /
C DATA LOG10(3),LOG10(4) / O102373, O147770 /
C
C MACHINE CONSTANTS FOR THE UNIVAC 1100 SERIES. FTN COMPILER
C
C DATA SMALL(1),SMALL(2) / O000040000000, O000000000000 /
C DATA LARGE(1),LARGE(2) / O377777777777, O777777777777 /
C DATA RIGHT(1),RIGHT(2) / O170540000000, O000000000000 /
C DATA DIVER(1),DIVER(2) / O170640000000, O000000000000 /
C DATA LOG10(1),LOG10(2) / O177746420232, O411757177572 /
C
C MACHINE CONSTANTS FOR THE UNIVAC 1100 SERIES. FOR COMPILER
C
C DATA SMALL(1), SMALL(2) / O000040000000, O000000000000 /
C DATA LARGE(1), LARGE(2) / O377777777777, O777777777777 /
C DATA RIGHT(1), RIGHT(2) / O170540000000, O000000000000 /
C DATA DIVER(1), DIVER(2) / O170640000000, O000000000000 /
C DATA LOG10(1), LOG10(2) / O177746420232, O411757177572/
C
C
C MACHINE CONSTANTS FOR VAX 11/780
C (EXPRESSED IN INTEGER AND HEXADECIMAL)
C ***THE HEX FORMAT BELOW MAY NOT BE SUITABLE FOR UNIX SYSYEMS***
C *** THE INTEGER FORMAT SHOULD BE OK FOR UNIX SYSTEMS***
C
C DATA SMALL(1), SMALL(2) / 128, 0 /
C DATA LARGE(1), LARGE(2) / -32769, -1 /
C DATA RIGHT(1), RIGHT(2) / 9344, 0 /
C DATA DIVER(1), DIVER(2) / 9472, 0 /
C DATA LOG10(1), LOG10(2) / 546979738, -805665541 /
C
C DATA SMALL(1), SMALL(2) / Z00000080, Z00000000 /
C DATA LARGE(1), LARGE(2) / ZFFFF7FFF, ZFFFFFFFF /
C DATA RIGHT(1), RIGHT(2) / Z00002480, Z00000000 /
C DATA DIVER(1), DIVER(2) / Z00002500, Z00000000 /
C DATA LOG10(1), LOG10(2) / Z209A3F9A, ZCFFA84FB /
C
C***FIRST EXECUTABLE STATEMENT D1MACH
C
D1MACH = DMACH(I)
RETURN
C
END
REAL FUNCTION R1MACH(I)
C***BEGIN PROLOGUE R1MACH
C***DATE WRITTEN 790101 (YYMMDD)
C***REVISION DATE 820801 (YYMMDD)
C***CATEGORY NO. R1
C***KEYWORDS MACHINE CONSTANTS
C***AUTHOR FOX, P. A., (BELL LABS)
C HALL, A. D., (BELL LABS)
C SCHRYER, N. L., (BELL LABS)
C***PURPOSE RETURNS SINGLE PRECISION MACHINE DEPENDENT CONSTANTS
C***DESCRIPTION
C R1MACH CAN BE USED TO OBTAIN MACHINE-DEPENDENT PARAMETERS
C FOR THE LOCAL MACHINE ENVIRONMENT. IT IS A FUNCTION
C SUBROUTINE WITH ONE (INPUT) ARGUMENT, AND CAN BE CALLED
C AS FOLLOWS, FOR EXAMPLE
C
C A = R1MACH(I)
C
C WHERE I=1,...,5. THE (OUTPUT) VALUE OF A ABOVE IS
C DETERMINED BY THE (INPUT) VALUE OF I. THE RESULTS FOR
C VARIOUS VALUES OF I ARE DISCUSSED BELOW.
C
C SINGLE-PRECISION MACHINE CONSTANTS
C R1MACH(1) = B**(EMIN-1), THE SMALLEST POSITIVE MAGNITUDE.
C R1MACH(2) = B**EMAX*(1 - B**(-T)), THE LARGEST MAGNITUDE.
C R1MACH(3) = B**(-T), THE SMALLEST RELATIVE SPACING.
C R1MACH(4) = B**(1-T), THE LARGEST RELATIVE SPACING.
C R1MACH(5) = LOG10(B)
C***REFERENCES FOX, P.A., HALL, A.D., SCHRYER, N.L, *FRAMEWORK FOR
C A PORTABLE LIBRARY*, ACM TRANSACTIONS ON MATHE-
C MATICAL SOFTWARE, VOL. 4, NO. 2, JUNE 1978,
C PP. 177-188.
C***ROUTINES CALLED XERROR
C***END PROLOGUE R1MACH
C
INTEGER SMALL(2)
INTEGER LARGE(2)
INTEGER RIGHT(2)
INTEGER DIVER(2)
INTEGER LOG10(2)
C
REAL RMACH(5)
C
EQUIVALENCE (RMACH(1),SMALL(1))
EQUIVALENCE (RMACH(2),LARGE(1))
EQUIVALENCE (RMACH(3),RIGHT(1))
EQUIVALENCE (RMACH(4),DIVER(1))
EQUIVALENCE (RMACH(5),LOG10(1))
C
C MACHINE CONSTANTS FOR THE BURROUGHS 1700 SYSTEM.
C
C DATA RMACH(1) / Z400800000 /
C DATA RMACH(2) / Z5FFFFFFFF /
C DATA RMACH(3) / Z4E9800000 /
C DATA RMACH(4) / Z4EA800000 /
C DATA RMACH(5) / Z500E730E8 /
C
C MACHINE CONSTANTS FOR THE BURROUGHS 5700/6700/7700 SYSTEMS.
C
C DATA RMACH(1) / O1771000000000000 /
C DATA RMACH(2) / O0777777777777777 /
C DATA RMACH(3) / O1311000000000000 /
C DATA RMACH(4) / O1301000000000000 /
C DATA RMACH(5) / O1157163034761675 /
C
C MACHINE CONSTANTS FOR THE CDC 6000/7000 SERIES.
C
C DATA RMACH(1) / 00014000000000000000B /
C DATA RMACH(2) / 37767777777777777777B /
C DATA RMACH(3) / 16404000000000000000B /
C DATA RMACH(4) / 16414000000000000000B /
C DATA RMACH(5) / 17164642023241175720B /
C
C MACHINE CONSTANTS FOR THE CRAY 1
C
C DATA RMACH(1) / 200004000000000000000B /
C DATA RMACH(2) / 577777777777777777777B /
C DATA RMACH(3) / 377214000000000000000B /
C DATA RMACH(4) / 377224000000000000000B /
C DATA RMACH(5) / 377774642023241175720B /
C
C MACHINE CONSTANTS FOR THE DATA GENERAL ECLIPSE S/200
C
C NOTE - IT MAY BE APPROPRIATE TO INCLUDE THE FOLLOWING CARD -
C STATIC RMACH(5)
C
C DATA SMALL/20K,0/,LARGE/77777K,177777K/
C DATA RIGHT/35420K,0/,DIVER/36020K,0/
C DATA LOG10/40423K,42023K/
C
C MACHINE CONSTANTS FOR THE HARRIS 220
C
C DATA SMALL(1),SMALL(2) / [20000000, [00000201 /
C DATA LARGE(1),LARGE(2) / [37777777, [00000177 /
C DATA RIGHT(1),RIGHT(2) / [20000000, [00000352 /
C DATA DIVER(1),DIVER(2) / [20000000, [00000353 /
C DATA LOG10(1),LOG10(2) / [23210115, [00000377 /
C
C MACHINE CONSTANTS FOR THE HONEYWELL 600/6000 SERIES.
C
C DATA RMACH(1) / O402400000000 /
C DATA RMACH(2) / O376777777777 /
C DATA RMACH(3) / O714400000000 /
C DATA RMACH(4) / O716400000000 /
C DATA RMACH(5) / O776464202324 /
C
C MACHINE CONSTANTS FOR THE HP 2100
C
C 3 WORD DOUBLE PRECISION WITH FTN4
C
C DATA SMALL(1), SMALL(2) / 40000B, 1 /
C DATA LARGE(1), LARGE(2) / 77777B, 177776B /
C DATA RIGHT(1), RIGHT(2) / 40000B, 325B /
C DATA DIVER(1), DIVER(2) / 40000B, 327B /
C DATA LOG10(1), LOG10(2) / 46420B, 46777B /
C
C MACHINE CONSTANTS FOR THE HP 2100
C 4 WORD DOUBLE PRECISION WITH FTN4
C
C DATA SMALL(1), SMALL(2) / 40000B, 1 /
C DATA LARGE91), LARGE(2) / 77777B, 177776B /
C DATA RIGHT(1), RIGHT(2) / 40000B, 325B /
C DATA DIVER(1), DIVER(2) / 40000B, 327B /
C DATA LOG10(1), LOG10(2) / 46420B, 46777B /
C
C MACHINE CONSTANTS FOR THE IBM 360/370 SERIES,
C THE XEROX SIGMA 5/7/9, THE SEL SYSTEMS 85/86 AND
C THE PERKIN ELMER (INTERDATA) 7/32.
C
C DATA RMACH(1) / Z00100000 /
C DATA RMACH(2) / Z7FFFFFFF /
C DATA RMACH(3) / Z3B100000 /
C DATA RMACH(4) / Z3C100000 /
C DATA RMACH(5) / Z41134413 /
C
C MACHINE CONSTANTS FOR THE PDP-10 (KA OR KI PROCESSOR).
C
C DATA RMACH(1) / "000400000000 /
C DATA RMACH(2) / "377777777777 /
C DATA RMACH(3) / "146400000000 /
C DATA RMACH(4) / "147400000000 /
C DATA RMACH(5) / "177464202324 /
C
C MACHINE CONSTANTS FOR PDP-11 FORTRAN SUPPORTING
C 32-BIT INTEGERS (EXPRESSED IN INTEGER AND OCTAL).
C
C DATA SMALL(1) / 8388608 /
C DATA LARGE(1) / 2147483647 /
C DATA RIGHT(1) / 880803840 /
C DATA DIVER(1) / 889192448 /
C DATA LOG10(1) / 1067065499 /
C
C DATA RMACH(1) / O00040000000 /
C DATA RMACH(2) / O17777777777 /
C DATA RMACH(3) / O06440000000 /
C DATA RMACH(4) / O06500000000 /
C DATA RMACH(5) / O07746420233 /
C
C MACHINE CONSTANTS FOR PDP-11 FORTRAN SUPPORTING
C 16-BIT INTEGERS (EXPRESSED IN INTEGER AND OCTAL).
C
C DATA SMALL(1),SMALL(2) / 128, 0 /
C DATA LARGE(1),LARGE(2) / 32767, -1 /
C DATA RIGHT(1),RIGHT(2) / 13440, 0 /
C DATA DIVER(1),DIVER(2) / 13568, 0 /
C DATA LOG10(1),LOG10(2) / 16282, 8347 /
C
C DATA SMALL(1),SMALL(2) / O000200, O000000 /
C DATA LARGE(1),LARGE(2) / O077777, O177777 /
C DATA RIGHT(1),RIGHT(2) / O032200, O000000 /
C DATA DIVER(1),DIVER(2) / O032400, O000000 /
C DATA LOG10(1),LOG10(2) / O037632, O020233 /
C
C MACHINE CONSTANTS FOR THE UNIVAC 1100 SERIES.
C
C DATA RMACH(1) / O000400000000 /
C DATA RMACH(2) / O377777777777 /
C DATA RMACH(3) / O146400000000 /
C DATA RMACH(4) / O147400000000 /
C DATA RMACH(5) / O177464202324 /
C
C MACHINE CONSTANTS FOR THE VAX 11/780
C (EXPRESSED IN INTEGER AND HEXADECIMAL)
C ***THE HEX FORMAT BELOW MAY NOT BE SUITABLE FOR UNIX SYSTEMS***
C *** THE INTEGER FORMAT SHOULD BE OK FOR UNIX SYSTEMS***
C
C DATA SMALL(1) / 128 /
C DATA LARGE(1) / -32769 /
C DATA RIGHT(1) / 13440 /
C DATA DIVER(1) / 13568 /
C DATA LOG10(1) / 547045274 /
C
C DATA SMALL(1) / Z00000080 /
C DATA LARGE(1) / ZFFFF7FFF /
C DATA RIGHT(1) / Z00003480 /
C DATA DIVER(1) / Z00003500 /
C DATA LOG10(1) / Z209B3F9A /
C
C***FIRST EXECUTABLE STATEMENT R1MACH
C
R1MACH = RMACH(I)
RETURN
C
END
3 DL1BL3= BLISS VOL II P.294... ALL DATA
60 1 5.0
(1X,F6.3,F3.0,2F6.3)
125 1 300 402
164 1 156 258
-065 1 401 503
120 1 137 239
034 1 214 316
010 1 084 186
-079 1 137 239
-533 1 031 133 X
-611 1 227 329 X
076 1 092 194
319 1 427 228
465 1 487 288
313 1 487 288
193 1 487 288
129 1 589 390
209 1 571 372
284 1 446 247
173 1 570 371
264 1 472 273
132 1 476 277
324 1 696 321
367 1 729 354
321 1 509 134
375 1 559 184
345 1 679 304
341 1 583 208
327 1 742 367
256 1 781 406
256 1 739 364
214 1 865 490
501 1 723 223
318 1 940 440
317 1 903 403
349 1 910 410
402 1 684 184
374 1 904 404
340 1 887 387
598 1 593 093
444 1 640 140
543 1 512 012
559 1 832 156
705 1 817 141
585 1 881 205
588 1 848 172
532 1 857 181
593 1 808 132
593 1 829 153
559 1 872 196
611 1 805 129
561 1 904 228
561 1 1047 246
733 1 942 141
553 1 1194 393
593 1 1090 289
525 1 1000 199
580 1 1054 253
711 1 956 155
546 1 1101 300
632 1 1060 259
733 1 893 092
3 DL1BL3= BLISS VOL II P.294... ALL DATA
PROBLEM TITLE DL1BL3= BLISS VOL II P.294... ALL DATA
NUMBER OF OBSERVATION= 60
NUMBER OF PARAMETERS= 3
MINIMUM NUMBER OF PARAMETERS CONSIDERED= 1
PERCENTAGE DEVIATION FROM OPTIMALITY ALLOWED= 5.00
**BEST SUBSET LAV PROGRAM
IFAULT= 0
BEST RESULTS FOR LAV SUBSET OF SIZE= 1
**************SUM OF ABSOLUTE VALUES= 7.659
BETA( 2) .550
BEST RESULTS FOR LAV SUBSET OF SIZE= 2
**************SUM OF ABSOLUTE VALUES= 5.407
BETA( 2) .818
BETA( 3) -.782
BEST RESULTS FOR LAV SUBSET OF SIZE= 3
**************SUM OF ABSOLUTE VALUES= 3.877
BETA( 3) -1.116
BETA( 2) .628
BETA( 1) .235
ITERATION COUNT= 106
0